Electrohydrodynamic lubrication
cp
PATH
AMPL
short
= 1, := 100; /* grid 0 .. N */
/* half-points go 1 .. N */
param pi := 4 * atan(1);
param xa := -3;
param xf > xa, := 2;
param dx := (xf - xa) / N;
param alpha := 2.832;
param lambda := 6.057;
param k_init := 1.6;
param p_init {i in 1..N} := max(0, 1 - abs((xa +1 + i*dx)/2));
param w {i in 0..N} := (if i in {0,N} then 0.5 else 1);
var k := k_init;
var p {i in 1..N} >= 0, := p_init[i];
var h {j in 0..N} =
( xa + (j+.5)*dx)^2 + k + 1
+ 1/pi * (sum {i in 0..N} w[i] * (i-j-0.5)*dx
* log(abs(i-j-.5)*dx) *
((if i < N then p[i+1]) - (if i > 1 then p[i-1]) )
) ;
psum :
k complements
( 1 - dx*2/pi * sum {i in 1..N} w[i]*p[i] )
= 0;
reynolds {i in 1..N}:
p[i] >= 0 complements
( (lambda / dx) * (h[i]-h[i-1])
- (1/dx^2) * (
h[i]^3 * ((if i < N then p[i+1])-p[i]) /
exp(alpha*((if i < N then p[i+1])+p[i])*.5)
- h[i-1]^3 * (p[i]-(if i > 1 then p[i-1])) /
exp(alpha*(p[i]+(if i > 1 then p[i-1]))*.5)
)
)
>= 0;
]]>
solve;
printf "Pressure is:\n";
display p;